{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(regression)=\n",
    "# Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this chapter, you'll learn how to run linear regressions with code.\n",
    "\n",
    "If you're running this code (either by copying and pasting it, or by downloading it using the icons at the top of the page), you may need to the packages it uses by, for example, running `pip install packagename` on your computer's command line. (If you're not sure what a command line is, take a quick look at the basics of coding chapter.)\n",
    "\n",
    "Most of this chapter will rely on [statsmodels](https://www.statsmodels.org/stable/index.html) with some use of [**linearmodels**](https://bashtage.github.io/linearmodels/). Some of the material in this chapter follows [Grant McDermott](https://grantmcdermott.com/)'s excellent notes and the [Library of Statistical Translation](https://lost-stats.github.io/).\n",
    "\n",
    "### Notation and basic definitions\n",
    "\n",
    "Greek letters, like $\\beta$, are the truth and represent parameters. Modified Greek letters are an estimate of the truth, for example $\\hat{\\beta}$. Sometimes Greek letters will stand in for vectors of parameters. Most of the time, upper case Latin characters such as $X$ will represent random variables (which could have more than one dimension). Lower case letters from the Latin alphabet denote realised data, for instance $x$ (which again could be multi-dimensional).  Modified Latin alphabet letters denote computations performed on data, for instance $\\bar{x} = \\frac{1}{n} \\displaystyle\\sum_{i} x_i$ where $n$ is number of samples.\n",
    "\n",
    "Ordinary least squares (OLS) regression can be used to *estimate* the parameters of certain types of model, most typically models of the form\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 \\cdot x_1 + \\beta_2 \\cdot x_2\n",
    "$$\n",
    "\n",
    "This generic model says that the value of an outcome variable $y$ is a linear function of one or more input predictor variables $x_i$, where the $x_i$ could be transforms of original data. But the above equation is a platonic ideal, what we call a data generating process (DGP). OLS allows us to recover *estimates* of the parameters of the model , i.e. to find $\\hat{\\beta_i}$ and to enable us to write an estimated model:\n",
    "\n",
    "$$\n",
    "y = \\hat{\\beta_0} + \\hat{\\beta_1} \\cdot x_1 + \\hat{\\beta_2} \\cdot x_2 + \\epsilon\n",
    "$$\n",
    "\n",
    "This equation can also be expressed in matrix form as\n",
    "\n",
    "$$\n",
    "y = x'\\cdot \\hat{\\beta} + \\epsilon\n",
    "$$\n",
    "\n",
    "where $x' = (1, x_1, \\dots, x_{n})'$ and $\\hat{\\beta} = (\\hat{\\beta_0}, \\hat{\\beta_1}, \\dots, \\hat{\\beta_{n}})$.\n",
    "\n",
    "Given data $y_i$ stacked to make a vector $y$ and $x_{i}$ stacked to make a matrix $X$, this can be solved for the coefficients $\\hat{\\beta}$ according to\n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = \\left(X'X\\right)^{-1} X'y\n",
    "$$\n",
    "\n",
    "To be sure that the estimates of these parameters are the *best linear unbiased estimate*, a few conditions need to hold: the Gauss-Markov conditions:\n",
    "\n",
    "1. $y$ is a linear function of the $\\beta_i$\n",
    "2. $y$ and the $x_i$ are randomly sampled from the population.\n",
    "3. There is no perfect multi-collinearity of variables.\n",
    "4. $\\mathbb{E}(\\epsilon | x_1, \\dots, x_n) = 0$ (unconfoundedness)\n",
    "5. $\\text{Var}(\\epsilon | x_1, \\dots, x_n) = \\sigma^2$ (homoskedasticity)\n",
    "\n",
    "(1)-(4) also guarantee that OLS estimates are unbiased and $\\mathbb{E}(\\hat{\\beta}_i) = \\beta_i$.\n",
    "\n",
    "The classic linear model requires a 6th assumption; that $\\epsilon \\thicksim \\mathcal{N}(0, \\sigma^2)$.\n",
    "\n",
    "The interpretation of regression coefficients depends on what their units are to begin with, but you can always work it out by differentiating both sides of the model equation with respect to the $x_i$. For example, for the first model equation above\n",
    "\n",
    "$$\n",
    "\\frac{\\partial y}{\\partial x_i} = \\beta_i\n",
    "$$\n",
    "\n",
    "so we get the interpretation that $\\beta_i$ is the rate of change of y with respect to $x_i$. If $x_i$ and $y$ are in levels, this means that a unit increase in $x_i$ is associated with a $\\beta_i$ units increase in $y$. If the right-hand side of the model is $\\ln x_i$ then we get\n",
    "\n",
    "$$\n",
    "\\frac{\\partial y}{\\partial x_i} = \\beta_i \\frac{1}{x_i} \n",
    "$$\n",
    "\n",
    "with some abuse of notation, we can rewrite this as $\\partial y = \\beta_i \\partial x_i/x_i$, which says that a percent change in $x_i$ is associated with a $\\beta_i$ unit change in $y$. With a logged $y$ variable, it's a percent change in $x_i$ that is associated with a percent change in $y$, or $\\partial y/y = \\beta_i \\partial x_i/x_i$ (note that both sides of this equation are unitless in this case). Finally, another example that is important in practice is that of log differences, eg $y = \\beta_i (\\ln x_i - \\ln x_i')$. Again, we will abuse notation and say that this case may be represented as $\\partial y = \\beta_i (\\partial x_i/x_i - \\partial x_i'/x_i')$, i.e. the difference in two percentages, a *percentage point* change, in $x_i$ is associated with a $\\beta_i$ unit change in $y$.\n",
    "\n",
    "### Imports\n",
    "\n",
    "Let's import some of the packages we'll be using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import os\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 6)\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression basics\n",
    "\n",
    "There are two ways to run regressions in [**statsmodels**](https://www.statsmodels.org/stable/index.html); passing the data directly as objects, and using formulae. We'll see both but, just to get things started, let's use the formula API.\n",
    "\n",
    "We'll use the starwars dataset to run a regression of mass on height for star wars characters. This example borrows very heavily from notes by [Grant McDermott](https://grantmcdermott.com/). First, let's bring the dataset in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/data/starwars.csv\",\n",
    "    index_col=0,\n",
    ")\n",
    "# Look at first few rows\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, now let's do a regression using OLS and a formula that says our y-variable is mass and our regressor is height:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = smf.ols(\"mass ~ height\", data=df).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, where are the results!? They're stored in the object we created. To peek at them we need to call the summary function (and, for easy reading, I'll print it out too using `print`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we're seeing here are really several tables glued together. To just grab the coefficients in a tidy format, use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll have noticed that we got an intercept, even though we didn't specify one in the formula. **statsmodels** adds in an intercept by default because, most of the time, you will want one. To turn it off, add a `-1` at the end of the formula command, eg in this case you would call `smf.ols('mass ~ height -1', data=df).fit()`.\n",
    "\n",
    "The fit we got in the case with the intercept was pretty terrible; a low $R^2$ and both of our confidence intervals are large and contain zero. What's going on? If there's one adage in regression that's always worth paying attention to, it's *always plot your data*. Let's see what's going on here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sns.scatterplot(data=df, x=\"height\", y=\"mass\", s=200, ax=ax, legend=False, alpha=0.8)\n",
    "ax.annotate(\n",
    "    \"Jabba the Hutt\",\n",
    "    df.iloc[df[\"mass\"].idxmax()][[\"height\", \"mass\"]],\n",
    "    xytext=(0, -50),\n",
    "    textcoords=\"offset points\",\n",
    "    arrowprops=dict(\n",
    "        arrowstyle=\"fancy\",\n",
    "        color=\"k\",\n",
    "        connectionstyle=\"arc3,rad=0.3\",\n",
    "    ),\n",
    ")\n",
    "ax.set_ylim(0, None)\n",
    "ax.set_title(\"Always plot the data\", loc=\"left\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Oh dear, Jabba's been on the paddy frogs again, and he's a bit of different case. When we're estimating statistical relationships, we have all kinds of choices and should be wary about arbitrary decisions of what to include or exclude in case we fool ourselves about the generality of the relationship we are capturing. Let's say we knew that we weren't interested in Hutts though, but only in other species: in that case, it's fair enough to filter out Jabba and run the regression without this obvious outlier. We'll exclude any entry that contains the string 'Jabba' in the `name` column:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_outlier_free = smf.ols(\n",
    "    \"mass ~ height\", data=df[~df[\"name\"].str.contains(\"Jabba\")]\n",
    ").fit()\n",
    "print(results_outlier_free.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks a lot more healthy. Not only is the model explaining a *lot* more of the data, but the coefficients are now significant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Robust regression\n",
    "\n",
    "Filtering out data is one way to deal with outliers, but it's not the only one; an alternative is to use a regression technique that is robust to such outliers. **statsmodels** has a variety of robust linear models that you can read more about [here](https://www.statsmodels.org/stable/examples/notebooks/generated/robust_models_0.html). To demonstrate the general idea, we will run the regression again but using a robust method.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_robust = smf.rlm(\n",
    "    \"mass ~ height\", data=df, M=sm.robust.norms.TrimmedMean(0.5)\n",
    ").fit()\n",
    "print(results_robust.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many different 'M-estimators' available; in this case the TrimmedMean estimator gives a very similar result to the regression with the point excluded. We can visualise this, and, well, the results are not really very different in this case. Note that `abline_plot` just takes an intercept and coefficient from a fitted model and renders the line that they encode."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.scatter(df[\"height\"], df[\"mass\"])\n",
    "sm.graphics.abline_plot(model_results=results_robust, ax=ax, alpha=0.5, label=\"Robust\")\n",
    "sm.graphics.abline_plot(\n",
    "    model_results=results, ax=ax, color=\"red\", label=\"OLS\", alpha=0.5, ls=\"--\"\n",
    ")\n",
    "ax.legend()\n",
    "ax.set_xlabel(\"Height\")\n",
    "ax.set_ylabel(\"Mass\")\n",
    "ax.set_ylim(0, None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard errors\n",
    "\n",
    "You'll have seen that there's a column for the standard error of the estimates in the regression table and a message saying that the covariance type of these is 'nonrobust'. Let's say that, instead, we want to use Eicker-Huber-White robust standard errors, aka \"HC2\" standard errors. We can specify to use these up front standard errors up front in the fit method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(smf.ols(\"mass ~ height\", data=df).fit(cov_type=\"HC2\").summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, alternatively, we can go back to our existing results and recompute the results from those:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.get_robustcov_results(\"HC2\").summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several different types of standard errors available in **statsmodels**:\n",
    "\n",
    "- ‘HC0’, ‘HC1’, ‘HC2’, and ‘HC3’\n",
    "- ‘HAC’, for heteroskedasticity and autocorrelation consistent standard errors, for which you may want to also use some keyword arguments\n",
    "- 'hac-groupsum’, for Driscoll and Kraay heteroscedasticity and\n",
    "autocorrelation robust standard errors in panel data, again for which you may have to specify extra keyword arguments\n",
    "- 'hac-panel’, for heteroscedasticity and autocorrelation robust standard\n",
    "errors in panel data, again with keyword arguments; and\n",
    "- 'cluster' for clustered standard errors.\n",
    "\n",
    "You can find information on all of these [here](https://www.statsmodels.org/stable/generated/statsmodels.regression.linear_model.OLSResults.get_robustcov_results.html?highlight=get_robustcov_results#statsmodels.regression.linear_model.OLSResults.get_robustcov_results). For more on standard errors in python, [this is a good](http://www.vincentgregoire.com/standard-errors-in-python/) link.\n",
    "\n",
    "For now, let's look more closely at those last ones: clustered standard errors.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clustered standard errors\n",
    "\n",
    "Often, we know something about the structure of likely errors, namely that they occur in groups. In the below example we use one-way clusters to capture this effect in the errors.\n",
    "\n",
    "Note that in the below example, we grab a subset of the data for which a set of variables we're interested in are defined, otherwise the below example would execute with an error because of missing cluster-group values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf = df.dropna(subset=[\"homeworld\", \"mass\", \"height\", \"species\"])\n",
    "results_clus = smf.ols(\"mass ~ height\", data=xf).fit(\n",
    "    cov_type=\"cluster\", cov_kwds={\"groups\": xf[\"homeworld\"]}\n",
    ")\n",
    "print(results_clus.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add two-way clustering of standard errors using the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xf = df.dropna(subset=[\"homeworld\", \"mass\", \"height\", \"species\"])\n",
    "two_way_clusters = np.array(xf[[\"homeworld\", \"species\"]], dtype=str)\n",
    "results_clus = smf.ols(\"mass ~ height\", data=xf).fit(\n",
    "    cov_type=\"cluster\", cov_kwds={\"groups\": two_way_clusters}\n",
    ")\n",
    "print(results_clus.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you would generally expect, the addition of clustering has increased the standard errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fixed effects and categorical variables\n",
    "\n",
    "Fixed effects are a way of allowing the intercept of a regression model to vary freely across individuals or groups. It is, for example, used to control for any individual-specific attributes that do not vary across time in panel data.\n",
    "\n",
    "Let's use the 'mtcars' dataset to demonstrate this. We'll read it in and set the datatypes of some of the columns at the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg = pd.read_csv(\n",
    "    \"https://raw.githubusercontent.com/LOST-STATS/lost-stats.github.io/source/Data/mtcars.csv\",\n",
    "    dtype={\"model\": str, \"mpg\": float, \"hp\": float, \"disp\": float, \"cyl\": \"category\"},\n",
    ")\n",
    "mpg.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have our data in we want to regress mpg (miles per gallon) on hp (horsepower) with fixed effects for cyl (cylinders). Now we *could* just pop in a formula like this `'mpg ~ hp + cyl'` because we took the trouble to declare that `cyl` was of datatype category when reading it in from the csv file. This means that statsmodels will treat it as a category and use it as a fixed effect by default.\n",
    "\n",
    "But when I read that formula I get nervous that `cyl` might not have been processed correctly (ie it could have been read in as a float, which is what it looks like) and it might just be treated as a float (aka a continuous variable) in the regression. Which is not what we want at all. So, to be safe, and make our intentions explicit (even when the data is of type 'category'), it's best to use the syntax `C(cyl)` to ask for a fixed effect.\n",
    "\n",
    "Here's a regression which does that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_fe = smf.ols(\"mpg ~ hp + C(cyl)\", data=mpg).fit()\n",
    "print(results_fe.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see here that two of the three possible values of `cyl`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"cyl\"].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "have been added as fixed effects regressors. The way that `+C(cyl)` has been added makes it so that the coefficients given are relative to the coefficient for the intercept. We can turn the intercept off to get a coefficient per unique `cyl` value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(smf.ols(\"mpg ~ hp + C(cyl) -1\", data=mpg).fit().summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When there is an intercept, the coefficients of fixed effect variables can be interpreted as being the average of $y$ for that class *compared* to the excluded classes holding all other categories and variables fixed. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### High dimensional fixed effects, aka absorbing regression\n",
    "\n",
    "Sometimes, you just have a LOT of fixed effects (and perhaps you don't particularly care about them individually). A common example is having a large number of firms as part of a panel. Fortunately, there are ways to make regression with high dimensional fixed effects be both fast and concise. (In Stata, this is provided by the `reghdfe` package.) Here, we will use the [**linearmodels**](https://bashtage.github.io/linearmodels/index.html) package, which is built on top of **statsmodels**.\n",
    "\n",
    "Let's say we have a regression of the form\n",
    "\n",
    "$$\n",
    "y_i = x_i\\cdot \\beta + z_i\\cdot \\gamma +\\epsilon_i\n",
    "$$\n",
    "\n",
    "where $y_i$ are observations indexed by $i$, $x_i$ are vectors of exogenous variables we care about the coefficients ($\\beta$), $z_i$ are vectors of fixed effects we don't care too much about the coefficients (\\gamma) for, and the $\\epsilon_i$ are errors. Then we can use an *absorbing regression* to solve for the $\\beta$ while ignoring the $\\gamma$.\n",
    "\n",
    "Here's an example using simulated data on workers taken from the **linearmodels** docs. Let's simulate some data first, with two fixed effects (state and firm) alongside the two exogenous variables we're interested in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import default_rng\n",
    "\n",
    "rng = default_rng()  # Random number generator\n",
    "\n",
    "# Create synthetic input data\n",
    "nobs = 1_000_000  # No. observations\n",
    "state_id = rng.integers(50, size=nobs)  # State identifier\n",
    "firm_id = rng.integers(nobs // 5, size=nobs)  # Firm identifier (mean of 5 workers/firm)\n",
    "x = rng.standard_normal((nobs, 2))  # Exogenous variables\n",
    "sim = pd.DataFrame(\n",
    "    {\n",
    "        \"state_id\": pd.Categorical(state_id),\n",
    "        \"firm_id\": pd.Categorical(firm_id),\n",
    "        \"exog_0\": x[:, 0],\n",
    "        \"exog_1\": x[:, 1],\n",
    "    }\n",
    ")\n",
    "\n",
    "# Create synthetic relationship\n",
    "beta = [1, 3]  # coefficients of interest\n",
    "state_effects = rng.standard_normal(state_id.max() + 1)\n",
    "state_effects = state_effects[state_id]  # Generate state fixed effects\n",
    "firm_effects = rng.standard_normal(firm_id.max() + 1)\n",
    "firm_effects = firm_effects[firm_id]  # Generate firm fixed effects\n",
    "eps = rng.standard_normal(nobs)  # Generate errors\n",
    "# Generate endogeneous outcome variable\n",
    "sim[\"y\"] = (\n",
    "    sim[\"exog_0\"] * beta[0]\n",
    "    + sim[\"exog_1\"] * beta[1]\n",
    "    + firm_effects\n",
    "    + state_effects\n",
    "    + eps\n",
    ")\n",
    "sim.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we pass this to **linearmodels** and with the `state_id` and `firm_id` variables entered via the `absorb` keyword argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.iv.absorbing import AbsorbingLS\n",
    "\n",
    "mod = AbsorbingLS(\n",
    "    sim[\"y\"], sim[[\"exog_0\", \"exog_1\"]], absorb=sim[[\"state_id\", \"firm_id\"]]\n",
    ")\n",
    "print(mod.fit())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, from our 1,000,000 observations, we have roughly 200,000 fixed effects that have been scooped up and packed away, leaving us with just the coefficients, $\\beta$, on the exogenous variables of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations of regressors\n",
    "\n",
    "This chapter is showcasing *linear* regression. What that means is that the model is linear in the regressors: but it doesn't mean that those regressors can't be some kind of (potentially non-linear) transform of the original features $x_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Logs and arcsinh\n",
    "\n",
    "You have two options for adding in logs: do them before, or do them in the formula. Doing them before just makes use of standard dataframe operations to declare a new column:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"lnhp\"] = np.log(mpg[\"hp\"])\n",
    "print(smf.ols(\"mpg ~ lnhp\", data=mpg).fit().summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can specify the log directly in the formula:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_ln = smf.ols(\"mpg ~ np.log(hp)\", data=mpg).fit()\n",
    "print(results_ln.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, the first method will work for `arcsinh(x)` and `log(x+1)`, but you can also pass both of these into the formula directly too. (For more on the pros and cons of arcsinh, see {cite}`bellemare2020elasticities`.) Here it is with arcsinh:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(smf.ols(\"mpg ~ np.arcsinh(hp)\", data=mpg).fit().summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interaction terms and powers\n",
    "\n",
    "This chapter is showcasing *linear* regression. What that means is that the model is linear in the regressors: but it doesn't mean that those regressors can't be some kind of non-linear transform of the original features $x_i$. Two of the most common transformations that you might want to use are *interaction terms* and *polynomial terms*. An example of an interaction term would be\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 x_1 \\cdot x_2\n",
    "$$\n",
    "\n",
    "while an example of a polynomial term would be\n",
    "\n",
    "$$\n",
    "y = \\beta_0 + \\beta_1 x_1^2\n",
    "$$\n",
    "\n",
    "i.e. the last term enters only after it is multiplied by itself.\n",
    "\n",
    "One note of warning: the interpretation of the effect of a variable is no longer as simple as was set out at the start of this chapter. To work out *what* the new interpretation is, the procedure is the same though: just take the derivative. In the case of the interaction model above, the effect of a unit change in $x_1$ on $y$ is now going to be a function of $x_2$. In the case of the polynomial model above, the effect of a unit change in $x_1$ on $y$ will be $2\\beta_1 \\cdot x_1$. For more on interaction terms, see {cite}`balli2013interaction`.\n",
    "\n",
    "Alright, with all of that preamble out of the way, let's see how we actual do some of this! Let's try including a linear and squared term in the regression of `mpg` on `hp` making use of the numpy power function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_poly = smf.ols(\"mpg ~ hp + np.power(hp, 2)\", data=mpg).fit()\n",
    "print(res_poly.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's include the original term in hp, a term in disp, and the interaction between them, which is represented by hp:disp in the table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_inter = smf.ols(\"mpg ~ hp * disp\", data=mpg).fit()\n",
    "print(res_inter.summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the unusual case that you want *only* the interaction term, you write it as it appears in the table above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(smf.ols(\"mpg ~ hp : disp\", data=mpg).fit().summary().tables[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The formula API explained\n",
    "\n",
    "As you will have seen `~` separates the left- and right-hand sides of the regression. `+` computes a set union, which will also be familiar from the examples above (ie it inludes two terms as long as they are distinct). `-` computes a set difference; it adds the set of terms to the left of it while removing any that appear on the right of it. As we've seen, `a*b` is a short-hand for `a + b + a:b`, with the last term representing the interaction. `/` is short hand for `a + a:b`, which is useful if, for example `b` is nested within `a`, so it doesn't make sense to control for `b` on its own. Actually, the `:` character can interact multiple terms so that `(a + b):(d + c)` is the same as `a:c + a:d + b:c + b:d`. `C(a)` tells statsmodels to treat `a` as a categorical variable that will be included as a fixed effect. Finally, as we saw above with powers, you can also pass in vectorised functions, such as `np.log` and `np.power`, directly into the formulae.\n",
    "\n",
    "One gotcha with the formula API is ensuring that you have sensible variable names in your dataframe, i.e. ones that do *not* include whitespace or, to take a really pathological example, have the name 'a + b' for one of the columns that you want to regress on. You can dodge this kind of problem by passing in the variable name as, for example, `Q(\"a + b\")` to be clear that the *column name* is anything within the `Q(\"...\")`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiple regression models\n",
    "\n",
    "As is so often the case, you're likely to want to run more than one model at once with different specifications. Although there is a base version of this in **statsmodels**, called `summary_col`, which you can find an example of [here](http://aeturrell.com//2018/05/05/running-many-regressions-alongside-pandas/), instead we'll be using the [**stargazer**](https://github.com/mwburke/stargazer) package to assemble the regressions together in a table.\n",
    "\n",
    "In the above examples, we've collected a few different regression results. Let's put them together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stargazer.stargazer import Stargazer\n",
    "\n",
    "\n",
    "stargazer_tab = Stargazer([results_ln, res_poly, res_inter])\n",
    "stargazer_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are lots of customisation options, including ones that add a title, rename variables, add notes, and so on. What is most useful is that as well as the HTML friendly output that you can see above, the package also exports to latex:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(stargazer_tab.render_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And of course this can be written to a file using `open('regression.tex', 'w').write(stargazer.render_latex())` where you can get your main latex compilation to scoop it up and use it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying regressions without formulae, using the array API\n",
    "\n",
    "As noted, there are two ways to run regressions in [**statsmodels**](https://www.statsmodels.org/stable/index.html); passing the data directly as objects, and using formulae. We've seen the formula API, now let's see how to specify regressions using arrays with the format `sm.OLS(y, X)`.\n",
    "\n",
    "We will first need to take the data out of the **pandas** dataframe and put it into a couple of arrays. When we're not using the formula API, the default is to treat the array X as the design matrix for the regression-so, if it doesn't have a column of constants in, there will be no intercept in the regression. Therefore, we need to add a constant vector to the matrix `X` if we *do* want an intercept. Use `sm.add_constant(X)` for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(xf[\"height\"])\n",
    "y = np.array(xf[\"mass\"])\n",
    "X = sm.add_constant(X)\n",
    "results = sm.OLS(y, X).fit()\n",
    "print(results.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "This approach seems a lot less convenient, not to mention less clear, so you may be wondering when it is useful. It's useful when you want to do many regressions in a systematic way or when you don't know what the columns of a dataset will be called ahead of time. It can actually be a little bit simpler to specify for more complex regressions too."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fixed effects in the array API\n",
    "\n",
    "If you're using the formula API, it's easy to turn a regressor `x` into a fixed effect by putting `C(x)` into the model formula, as you'll see in the next section.\n",
    "\n",
    "For the array API, things are not that simple and you need to use dummy variables. Let's say we have some data like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import Generator, PCG64\n",
    "\n",
    "# Set seed for random numbers\n",
    "seed_for_prng = 78557\n",
    "prng = Generator(PCG64(seed_for_prng))\n",
    "\n",
    "no_obs = 200\n",
    "X = pd.DataFrame(prng.normal(size=no_obs))\n",
    "X[1] = prng.choice([\"a\", \"b\"], size=no_obs)\n",
    "# Get this a numpy array\n",
    "X = X.values\n",
    "# Create the y data, adding in a bit of noise\n",
    "y = X[:, 0] * 2 + 0.5 + prng.normal(scale=0.1, size=no_obs)\n",
    "y = [el_y + 1.5 if el_x == \"a\" else el_y + 3.4 for el_y, el_x in zip(y, X[:, 1])]\n",
    "X[:5, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first feature (column) is of numbers and it's clear how we include it. The second, however, is a grouping that we'd like to include as a fixed effect. But if we just throw this matrix into `sm.OLS(y, X)`, we're going to get trouble because **statsmodels** isn't sure what to do with a vector of strings. So, instead, we need to create some dummy variables out of our second column of data\n",
    "\n",
    "Astonishingly, there are several popular ways to create dummy variables in Python: **scikit-learn**'s `OneHotEncoder` and **pandas**' `get_dummies` being my favourites. Let's use the latter here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.get_dummies(X[:, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just need to pop this into our matrix $X$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.column_stack([X[:, 0], pd.get_dummies(X[:, 1])])\n",
    "X = np.array(X, dtype=float)\n",
    "X[:5, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, so now we're ready to do our regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sm.OLS(y, X).fit().summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perhaps you can see why I generally prefer the formula API..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instrumental variables\n",
    "\n",
    "Rather than use **statsmodels** for IV, we'll use the [**linearmodels**](https://bashtage.github.io/linearmodels/doc/index.html) package, which has very clean documentation (indeed, this sub-section is indebted to that documentation).\n",
    "\n",
    "\n",
    "Recall that a good instrumental variable $z$ has zero covariance with the error from the regression (which is untestable) and non-zero covariance with the variable of interest (which is).\n",
    "\n",
    "\n",
    "Recall that in IV regression, we have a model of the form\n",
    "\n",
    "$$\n",
    "\\begin{split}y_i    & = x_{1i}\\hat{\\beta_1} + x_{2i}\\hat{\\beta_2} + \\epsilon_i \\\\\n",
    "x_{2i} & = z_{1i}\\hat{\\delta} + z_{2i}\\hat{\\gamma} + \\nu_i\\end{split}\n",
    "$$\n",
    "\n",
    "where $x_{1i}$ is a set of $k_1$ exogenous regressors and $x_{2i}$ is a set of $k_2$ endogenous regressors such that $\\text{Cov}(x_{2i}, \\epsilon_i)\\neq 0$. This is a problem for the usual OLS assumptions (the right-hand side should be exogenous). \n",
    "\n",
    "\n",
    "To get around this, in 2-stage least squares IV, we first regress $x_{2i}$ on instruments that explain $x_{2i}$ *but not* $y_i$, and then regress $y_i$ only on the predicted/estimated left-hand side from the first regression, ie on $\\hat{x_{2i}}$. There are other estimators than IV2SLS, but I think that one has the most intuitive explanation of what's going.\n",
    "\n",
    "As well as a 2-stage least squares estimator called `IV2SLS`, **linearmodels** has a Limited Information Maximum Likelihood (LIML) estimator `IVLIML`, a Generalized Method of Moments (GMM) estimator `IVGMM`, and a Generalized Method of Moments using the Continuously Updating Estimator (CUE) `IVGMMCUE`.\n",
    "\n",
    "Just as with OLS via **statsmodels**, there's an option to use an array API for the **linearmodels** IV methods.\n",
    "\n",
    "It's always easiest to see an example, so let's estimate what might cause (realised) cigarette demand for the 48 continental US states in 1995 with `IV2SLS`. First we need to import the estimator, `IV2SLS`, and the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linearmodels.iv import IV2SLS\n",
    "\n",
    "df = pd.read_csv(\n",
    "    \"https://vincentarelbundock.github.io/Rdatasets/csv/AER/CigarettesSW.csv\",\n",
    "    dtype={\"state\": \"category\", \"year\": \"category\"},\n",
    ").assign(\n",
    "    rprice=lambda x: x[\"price\"] / x[\"cpi\"],\n",
    "    rincome=lambda x: x[\"income\"] / x[\"population\"] / x[\"cpi\"],\n",
    ")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll specify the model. It's going to be in the form `dep ~ exog + [endog ~ instruments]`, where endog will be regressed on instruments and dep will be regressed on both exog and the predicted values of endog.\n",
    "\n",
    "In this case, the model will be\n",
    "\n",
    "$$\n",
    "\\text{Price}_i = \\hat{\\pi_0} + \\hat{\\pi_1} \\text{SalesTax}_i + v_i \n",
    "$$\n",
    "\n",
    "in the first stage regression and\n",
    "\n",
    "$$\n",
    "\\text{Packs}_i = \\hat{\\beta_0} + \\hat{\\beta_2}\\widehat{\\text{Price}_i} + \\hat{\\beta_1} \\text{RealIncome}_i + u_i\n",
    "$$\n",
    "\n",
    "in the second stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_iv2sls = IV2SLS.from_formula(\n",
    "    \"np.log(packs) ~ 1 + np.log(rincome) + C(year) + C(state) + [np.log(rprice) ~ taxs]\",\n",
    "    df,\n",
    ").fit(cov_type=\"clustered\", clusters=df[\"year\"])\n",
    "print(results_iv2sls.summary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We sort of skipped a step here and did everything all in one go. If we *did* want to know how our first stage regression went, we can just pass a formula to `IV2SLS` without the part in square brackets, `[...]`, and it will run regular OLS. \n",
    "\n",
    "But, in this case, there's an easier way: we can print out a set of handy 1st stage statistics from running the full model.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_iv2sls.first_stage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are more tests and checks available. For example, Wooldridge’s regression test of exogeneity uses regression residuals from the endogenous variables regressed on the exogenous variables and the instrument to test for endogenity and is available to run on fitted model results. Let's check that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_iv2sls.wooldridge_regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the IV results against (naive) OLS. First, run the OLS equivalent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_cig_ols = IV2SLS.from_formula(\n",
    "    \"np.log(packs) ~ 1 + np.log(rincome) + C(year) + C(state) + np.log(rprice)\", df\n",
    ").fit(cov_type=\"clustered\", clusters=df[\"year\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now select these two models to compare:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "from linearmodels.iv.results import compare\n",
    "\n",
    "res = OrderedDict()\n",
    "res[\"OLS\"] = res_cig_ols\n",
    "res[\"2SLS\"] = results_iv2sls\n",
    "\n",
    "print(compare(res))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we take into account the fact that the real price is endogeneous to (realised) demand, we find that its coefficient is more negative; i.e. an increase in the real price of cigarettes creates a bigger fall in number of packs bought."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logit, probit, and generalised linear models\n",
    "\n",
    "### Logit\n",
    "\n",
    "A logistical regression, aka a logit, is a statistical method for a best-fit line between a regressors $X$ and an outcome varibale $y$ that takes on values in $(0, 1)$.\n",
    "\n",
    "The function that we're assuming links the regressors and the outcome has a few different names but the most common is the sigmoid function or the logistic function. The data generating process is assumed to be\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mathbb{P}(Y=1\\mid X) = \\frac{1}{1 + e^{-X'\\beta}}}\n",
    "$$\n",
    "\n",
    "we can also write this as $\\ln\\left(\\frac{p}{p-1}\\right) = \\beta_0 + \\sum_i \\beta_i x_i$ to get a 'log-odds' relationship. The coefficients from a logit model do not have the same interpration as in an OLS estimation, and you can see this from the fact that $\\partial y/\\partial x_i \\neq \\beta_i$ for logit. Of course, you can work out what the partial derivative is for yourself but most packages offer a convenient way to quickly recover the marginal effects. \n",
    "\n",
    "Logit models are available in **scikit-learn** and **statsmodels** but bear in mind that the **scikit-learn** logit model is, ermm, extremely courageous in that regularisation is applied by default. If you don't know what that means, don't worry, but it's probably best to stick with **statsmodels** as we will do in this example.\n",
    "\n",
    "We will predict a target `GRADE`, representing whether a grade improved or not, based on some regressors including participation in a programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data from Spector and Mazzeo (1980)\n",
    "df = sm.datasets.spector.load_pandas().data\n",
    "# Look at info on data\n",
    "print(sm.datasets.spector.NOTE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_logit = smf.logit(\"GRADE ~ GPA + TUCE + PSI\", data=df).fit()\n",
    "print(res_logit.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, did participation (`PSI`) help increase a grade? Yes. But we need to check the marginal effect to say exactly how much. We'll use `get_margeff` to do this, we'd like the $dy/dx$ effect, and we'll take it at the mean of each regressor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marg_effect = res_logit.get_margeff(at=\"mean\", method=\"dydx\")\n",
    "marg_effect.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So participation gives almost half a grade increase."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probit\n",
    "\n",
    "Probit is very similar to logit: it's a statistical method for a best-fit line between regressors $X$ and an outcome varibale $y$ that takes on values in $(0, 1)$. And, just like with logit, the function that we're assuming links the regressors and the outcome has a few different names! \n",
    "\n",
    "The data generating process is assumed to be\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mathbb{P}(Y=1\\mid X)=\\Phi (X^{T}\\beta )}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\Phi (x)={\\frac {1}{\\sqrt {2\\pi }}}\\int _{-\\infty }^{x}e^{-{\\frac {y^{2}}{2}}}dy.}\n",
    "$$\n",
    "\n",
    "is the cumulative standard normal (aka Gaussian) distribution. The coefficients from a probit model do not have the same interpration as in an OLS estimation, and you can see this from the fact that $\\partial y/\\partial x_i \\neq \\beta_i$ for probit. And, just as with logit, although you can derive the marginal effects, most packages offer a convenient way to quickly recover them. \n",
    "\n",
    "We can re-use our previous example of predicting a target `GRADE`, representing whether a grade improved or not, based on some regressors including participation (PSI) in a programme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_probit = smf.probit(\"GRADE ~ GPA + TUCE + PSI\", data=df).fit()\n",
    "print(res_probit.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_marg_effect = res_probit.get_margeff(at=\"mean\", method=\"dydx\")\n",
    "p_marg_effect.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's no coincidence that we find very similar results here because the two functions we're using don't actually look all that different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "support = np.linspace(-6, 6, 1000)\n",
    "ax.plot(support, st.logistic.cdf(support), \"r-\", ls=\"--\", label=\"Logistic\")\n",
    "ax.plot(support, st.norm.cdf(support), label=\"Probit\")\n",
    "ax.legend()\n",
    "ax.set_ylim(0, None)\n",
    "ax.set_ylim(0, None)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What difference there is, is that logistic regression puts more weight into the tails of the distribution. Arguably, logit is easier to interpret too. With logistic regression, a one unit change in $x_i$ is associated with a $\\beta_i$ change in the log odds of a 1 outcome or, alternatively, an $e^{\\beta_i}$-fold change in the odds, all else being equal. With a probit, this is a change of $\\beta_i z$ for $z$ a normalised variable that you'd have to convert into a predicted probability using the normal CDF."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generalised linear models\n",
    "\n",
    "Logit and probit (and OLS for that matter) as special cases of a class of models such that $g$ is a 'link' function connects a function of regressors to the output, and $\\mu$ is the mean of a conditional response distribution at a given point in the space of regressors. When $g(\\mu) = X'\\beta$, we just get regular OLS. When it's logit, we have\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\mu= \\mathbb{E}(Y\\mid X=x) =g^{-1}(X'\\beta)= \\frac{1}{1 + e^{-X'\\beta}}.}\n",
    "$$\n",
    "\n",
    "But as well as the ones we've seen, there are many possible link functions one can use via the catch-all `glm` function. These come in different 'families' of distributions, with the default for the binomial family being logit. So, running `smf.glm('GRADE ~ GPA + TUCE + PSI', data=df, family=sm.families.Binomial()).fit()` will produce exactly the same as we got both using the `logit` function. For more on the families of distributions and possible link functions, see the [relevant part](https://www.statsmodels.org/stable/glm.html#) of the **statsmodels** documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear probability model\n",
    "\n",
    "When $y$ takes values in $\\{0, 1\\}$ but the model looks like\n",
    "\n",
    "$$\n",
    "y = x' \\cdot \\beta\n",
    "$$\n",
    "\n",
    "and is estimated by OLS then you have a linear probability model. In this case, the interpretion of a unit change in $x_i$ is that it induces a $\\beta_i$ *change in probability* of $y$. Note that homoskedasticity does not hold for the linear probability model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Violations of the classical linear model (CLM)\n",
    "\n",
    "### Heteroskedasticity\n",
    "\n",
    "If an estimated model is homoskedastic then its random variables have equal (finite) variance. This is also known as homogeneity of variance. Another way of putting it is that, for all *observations* $i$ in an estimated model $y_i = X_i\\hat{\\beta} + \\epsilon_i$ then\n",
    "\n",
    "$$\n",
    "\\mathbb{E}(\\epsilon_i \\epsilon_i) = \\sigma^2\n",
    "$$\n",
    "\n",
    "When this relationship does not hold, an estimated model is said to be heteroskedastic.\n",
    "\n",
    "To test for heteroskedasticity, you can use **statsmodels**' versions of the [Breusch-Pagan](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.het_breuschpagan.html#statsmodels.stats.diagnostic.het_breuschpagan) or [White](https://www.statsmodels.org/stable/generated/statsmodels.stats.diagnostic.het_white.html#statsmodels.stats.diagnostic.het_white) tests with the null hypothesis that the estimated model is homoskedastic. If the null hypothesis is rejected, then standard errors, t-statistics, and F-statistics are invalidated. In this case, you will need HAC (heteroskedasticity and auto-correlation consistent) standard errors, t- and F-statistics.\n",
    "\n",
    "To obtain HAC standard errors from existing regression results in a variable `results`, you can use (for 1 lag):\n",
    "\n",
    "```python\n",
    "results.get_robustcov_results('HAC', maxlags=1).summary()\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantile regression\n",
    "\n",
    "Quantile regression estimates the conditional quantiles of a response variable. In some cases, it can be more robust to outliers and, in the case of the $q=0.5$ quantile it is equivalent LAD (Least Absolute Deviation) regression. Let's look at an example of quantile regression in action, lifted direct from the **statsmodels** [documentation](https://www.statsmodels.org/dev/examples/notebooks/generated/quantile_regression.html) and based on a Journal of Economic Perspectives paper by Koenker and Hallock."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = sm.datasets.engel.load_pandas().data\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have here are two sets of related data. Let's perform several quantile regressions from 0.1 to 0.9 in steps of 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mod = smf.quantreg(\"foodexp ~ income\", df)\n",
    "quantiles = np.arange(0.1, 1.0, 0.1)\n",
    "q_results = [mod.fit(q=x) for x in quantiles]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $q=0.5$ entry will be at the `4` index; let's take a look at it: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(q_results[4].summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the results for all of the regressions *and* let's add in OLS for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ols_res = smf.ols(\"foodexp ~ income\", df).fit()\n",
    "\n",
    "get_y = lambda a, b: a + b * x\n",
    "x = np.arange(df.income.min(), df.income.max(), 50)\n",
    "# Just to make the plot clearer\n",
    "x_max = 3000\n",
    "x = x[x < x_max]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "df.plot.scatter(\n",
    "    ax=ax, x=\"income\", y=\"foodexp\", alpha=0.7, s=10, zorder=2, edgecolor=None\n",
    ")\n",
    "for i, res in enumerate(q_results):\n",
    "    y = get_y(res.params[\"Intercept\"], res.params[\"income\"])\n",
    "    ax.plot(x, y, color=\"grey\", lw=0.5, zorder=0, linestyle=(0, (5, 10)))\n",
    "    ax.annotate(f\"$q={quantiles[i]:1.1f}$\", xy=(x.max(), y.max()))\n",
    "y = get_y(ols_res.params[\"Intercept\"], ols_res.params[\"income\"])\n",
    "ax.plot(x, y, color=\"red\", label=\"OLS\", zorder=0)\n",
    "ax.legend()\n",
    "ax.set_xlim(0, x_max)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chart shows very clearly how quantile regression differs from OLS. The line fitted by OLS is trying to be all things to all points whereas the line fitted by quantile regression is focused only on its quantile. You can also see how points far from the median (not all shown) may be having a large influence on the OLS line."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rolling and recursive regressions\n",
    "\n",
    "Rolling ordinary least squares applies OLS (ordinary least squares) across a fixed window of observations and then rolls (moves or slides) that window across the data set. They key parameter is `window`, which determines the number of observations used in each OLS regression. Recursive regression is equivalent to rolling regression but with a window that expands over time.\n",
    "\n",
    "Let's first create some synthetic data to perform estimation on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.regression.rolling import RollingOLS\n",
    "import statsmodels.api as sm\n",
    "from sklearn.datasets import make_regression\n",
    "\n",
    "X, y = make_regression(n_samples=200, n_features=2, random_state=0, noise=4.0, bias=0)\n",
    "df = pd.DataFrame(X).rename(columns={0: \"feature0\", 1: \"feature1\"})\n",
    "df[\"target\"] = y\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's fit the model using a formula and a `window` of 25 steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roll_reg = RollingOLS.from_formula(\n",
    "    \"target ~ feature0 + feature1 -1\", window=25, data=df\n",
    ")\n",
    "model = roll_reg.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that -1 in the formala suppresses the intercept. We can see the parameters using `model.params`. Here are the params for time steps between 20 and 30:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.params[20:30]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there aren't parameters for entries between 0 and 23 because our window is 25 steps wide. We can easily look at how any of the coefficients are changing over time. Here's an example for 'feature0'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = model.plot_recursive_coefficient(variables=[\"feature0\"])\n",
    "plt.xlabel(\"Time step\")\n",
    "plt.ylabel(\"Coefficient value\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A rolling regression with an *expanding* rather than *moving* window is effectively a recursive least squares model. We can do this instead using the `RecursiveLS` function from **statsmodels**. Let's fit this to the whole dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reg_rls = sm.RecursiveLS.from_formula(\"target ~ feature0 + feature1 -1\", df)\n",
    "model_rls = reg_rls.fit()\n",
    "print(model_rls.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But now we can look back at how the values of the coefficients changed over time too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = model_rls.plot_recursive_coefficient(\n",
    "    range(reg_rls.k_exog), legend_loc=\"upper right\"\n",
    ")\n",
    "ax_list = fig.axes\n",
    "for ax in ax_list:\n",
    "    ax.set_xlim(0, None)\n",
    "ax_list[-1].set_xlabel(\"Time step\")\n",
    "ax_list[0].set_title(\"Coefficient value\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Regression plots\n",
    "\n",
    "**statsmodels** has a number of built-in plotting methods to help you understand how well your regression is capturing the relationships you're looking for. Let's see a few examples of these using **statsmodels** built-in Statewide Crime data set:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crime_data = sm.datasets.statecrime.load_pandas()\n",
    "print(sm.datasets.statecrime.NOTE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's look at a Q-Q plot to get a sense of how the variables are distributed. This uses **scipy**'s stats module. The default distribution is normal but you can use any that **scipy** supports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "st.probplot(crime_data.data[\"murder\"], dist=\"norm\", plot=plt);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, this is not quite normal and there are some serious outliers in the tails.\n",
    "\n",
    "Let's run take a look at the unconditional relationship we're interested in: how murder depends on high school graduation. We'll use [**plotnine**](https://plotnine.readthedocs.io/en/stable/index.html)'s `geom_smooth` to do this but bear in mind it will only run a linear model of `'murder ~ hs_grad'` and ignore the other covariates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from plotnine import *\n",
    "\n",
    "(\n",
    "    ggplot(crime_data.data, aes(y=\"murder\", x=\"hs_grad\"))\n",
    "    + geom_point()\n",
    "    + geom_smooth(method=\"lm\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can take into account those other factors by using a partial regression plot that asks what does $\\mathbb{E}(y|X)$ look like as a function of $\\mathbb{E}(x_i|X)$? (Use `obs_labels=False` to remove data point labels.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"font.size\": 5}):\n",
    "    sm.graphics.plot_partregress(\n",
    "        endog=\"murder\",\n",
    "        exog_i=\"hs_grad\",\n",
    "        exog_others=[\"urban\", \"poverty\", \"single\"],\n",
    "        data=crime_data.data,\n",
    "        obs_labels=True,\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, the results of the regression are useful context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_crime = smf.ols(\n",
    "    \"murder ~ hs_grad + urban + poverty + single\", data=crime_data.data\n",
    ").fit()\n",
    "print(results_crime.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting the multicollinearity problems to one side, we see that the relationship shown in the partial regression plot is also implied by the coefficient on `hs_grad` in the regression table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also look at an in-depth summary of one exogenous regressor and its relationship to the outcome variable. Each of these types of regression diagnostic are available individually, or for all regressors at once, too. The first panel is the chart we did with **plotnine** rendered differently (and, one could argue, more informatively). Most of the plots below are self-explanatory except for the third one,  the CCPR (Component-Component plus Residual) plot. This provides a way to judge the effect of one regressor on the response variable by taking into account the effects of the other independent variables. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8, 6), dpi=150)\n",
    "\n",
    "sm.graphics.plot_regress_exog(results_crime, \"hs_grad\", fig=fig)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**statsmodels** can also produce influence plots of the 'externally studentised' residuals vs. the leverage of each observation as measured by the so-called hat matrix $X(X^{\\;\\prime}X)^{-1}X^{\\;\\prime}$ (because it puts the 'hat' on $y$). Externally studentised residuals are residuals that are scaled by their standard deviation. High leverage points could exert an undue influence over the regression line, but only if the predicted $y$ values of a regression that was fit with them excluded was quite different. In the example below, DC is having a big influence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({\"font.size\": 6}):\n",
    "    sm.graphics.influence_plot(results_crime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, it's nice to be able to see plots of our coefficients along with their standard errors. There isn't a built-in **statsmodels** option for this, but happily it's easy to extract the results of regressions in a sensible format. Using the `results` object from earlier, and excluding the intercept, we can get the coefficients from `results.params[1:]` and the associated errors from `results.bse[1:]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Put the results into a dataframe with Name, Coefficient, Error\n",
    "res_df = (\n",
    "    pd.concat([results_crime.params[1:], results_crime.bse[1:]], axis=1)\n",
    "    .reset_index()\n",
    "    .rename(columns={\"index\": \"Name\", 0: \"Coefficient\", 1: \"Error\"})\n",
    ")\n",
    "# Plot the coefficient values and their errors\n",
    "(\n",
    "    ggplot(res_df)\n",
    "    + geom_point(aes(\"Name\", \"Coefficient\"))\n",
    "    + geom_errorbar(aes(x=\"Name\", ymin=\"Coefficient-Error\", ymax=\"Coefficient+Error\"))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specification curve analysis\n",
    "\n",
    "When specifying a model, modellers have many options. These can be informed by field intelligence, priors, and even misguided attempts to find a significant result. Even with the best of intentions, research teams can reach entirely different conclusions using the same, or similar, data because of different choices made in preparing data or in modelling it.\n",
    "\n",
    "There’s formal evidence that researchers really do make different decisions; one study {cite}`silberzahn2018many` gave the same research question - whether soccer referees are more likely to give red cards to dark-skin-toned players than to light-skin-toned players - to 29 different teams. From the abstract of that paper:\n",
    "\n",
    "> Analytic approaches varied widely across the teams, and the estimated effect sizes ranged from 0.89 to 2.93 (Mdn = 1.31) in odds-ratio units. Twenty teams (69%) found a statistically significant positive effect, and 9 teams (31%) did not observe a significant relationship. Overall, the 29 different analyses used 21 unique combinations of covariates. Neither analysts’ prior beliefs about the effect of interest nor their level of expertise readily explained the variation in the outcomes of the analyses. Peer ratings of the quality of the analyses also did not account for the variability.\n",
    "\n",
    "So not only were different decisions made, there seems to be no clearly identifiable reason for them. There is usually scope for reasonable alternative model specifications when estimating coefficients, and those coefficients will vary with those specifications. \n",
    "\n",
    "Specification curve analysis {cite}`simonsohn2020specification` looks for a more exhaustive way of trying out alternative specifications. The three steps of specification curve analysis are:\n",
    "\n",
    "1. identifying the set of theoretically justified, statistically valid, and non-redundant analytic specifications;\n",
    "\n",
    "2. displaying alternative results graphically, allowing the identification of decisions producing different results; and\n",
    "\n",
    "3. conducting statistical tests to determine whether as a whole results are inconsistent with the null hypothesis.\n",
    "\n",
    "For a good example of specification curve analysis in action, see this recent Nature Human Behaviour paper {cite}`orben2019association` on the association between adolescent well-being and the use of digital technology.\n",
    "\n",
    "We'll use the [**specification curve analysis**](https://specification-curve.readthedocs.io/en/latest/readme.html) package to do the first two, which you can install with `pip install specification_curve` (full disclosure: I wrote this package). To demonstrate the full functionality, we'll create a second, alternative 'hp' that is a transformed version of the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mpg[\"hp_boxcox\"], _ = st.boxcox(mpg[\"hp\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create a specification curve. We need to specify the data, the different outcome variables we'd like to try, `y_endog`; the different possible versions of the main regressor of interest, `x_exog`; the possible controls, `controls`; any controls that should always be included, `always_include`; and any categorical variables to include class-by-class, `cat_expand`. Some of these accept lists of variables as well as single reggressors. The point estimates that have confidence intervals which include zero are coloured in grey, instead of blue. There is also an `exclu_grps` option to exclude certain combinations of regressors, and you can pass alternative estimators to fit, for example `fit(estimator=sm.Logit)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specification_curve import specification_curve as specy\n",
    "\n",
    "sc = specy.SpecificationCurve(\n",
    "    mpg,\n",
    "    y_endog=\"mpg\",\n",
    "    x_exog=[\"lnhp\", \"hp_boxcox\"],\n",
    "    controls=[\"drat\", \"qsec\", \"cyl\", \"gear\"],\n",
    "    always_include=[\"gear\"],\n",
    "    cat_expand=\"cyl\",\n",
    ")\n",
    "sc.fit()\n",
    "sc.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Review\n",
    "\n",
    "In this very short introduction to regression with code, you should have learned how to:\n",
    "\n",
    "- ✅ perform linear OLS regressions with code;\n",
    "- ✅ add fixed effects/categorical variables to regressions;\n",
    "- ✅ use different standard errors; \n",
    "- ✅ use models with transformed regressors;\n",
    "- ✅ use the formula or array APIs for **statsmodels** and **linearmodels**;\n",
    "- ✅ show the results from multiple models;\n",
    "- ✅ perform IV regressions;\n",
    "- ✅ perform GLM regressions; and\n",
    "- ✅ use plots as a way to interrogate regression results."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "codeforecon",
   "language": "python",
   "name": "codeforecon"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
